COUNTING ALL REGULAR OCTAHEDRONS IN {0,1,..., n} 3 

EUGEN J. IONASCU 

Abstract. In this paper we describe a procedure for calculating the number of regular octahedrons 
that have vertices with coordinates in the set {0, 1, ...,n}. As a result, we introduce a new sequence 
in The Online Encyclopedia of Integer Sequences (A178797) and list the first one hundred terms of it. 
We adapt the method appeared in [TT] which was used to find the number of regular tetrahedra with 
coordinates of their vertices in {0, 1, ..., n}. The idea of this calculation is based on the theoretical 
results obtained in [14] . A new fact proved here helps increasing the speed of all the programs used 
before. The procedure is put together in a series of commands written for Maple. 
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43 ■ 1. INTRODUCTION 

3- 

In this note we complete the work begun in the sequence of papers [2], [9]-|14j about equilat- 
eral triangles, regular tetrahedra, cubes, and regular octahedrons all with vertices having integer 
coordinates. Very often we will refer to this property by saying that the various objects are in Z 3 . 

£ I Strictly speaking these geometric objects are defined as being more than the set of their vertices 

that determines them, but for us here, these are just the vertices. So, for instance, an equilateral 
triangle is going to be a set of three points in Z 3 for which the Euclidean distances between every 
two of these points are the same. The main purpose of the paper is to take a close look at the 
regular octahedrons in Z . The simplest example of a regular octahedron with integer coordinates 
for its vertices is 

X 

S ; OCx := {[0, 1, 1], [1,0, 1], [1, 1, 0], [1, 1, 2], [1, 2, 1], [2, 1, 1]}, 

which can be obtained from the usual unit cube in R 3 , multiplying the vertices by a factor of two 
and then taking the coordinates of the centers of the new faces. It turns out that this procedure 
gives all such octahedrons as shown in |14j : 
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Figure 1 (a): OC\ octahedron Figure 1 (b): Regular octahedronvscube 

Theorem 1.1. Every regular octahedron in Z 3 is the dual of a cube that can be obtained (up 
to a translation with a vector with integer coordinates) by doubling a cube in Z 3 . 

Referring to Figure Q] (b), we showed that if the regular octahedron IJKLMN is in Z 3 , then so 
is the cube BB\C\IK\LOM and vice versa. This defines a one-to-one correspondence between the 
classes of cubes (invariant under integer translations) and the classes regular octahedra (invariant 
under integer translations) in Z 3 . In [I2j we determined a sequence of irreducible cubes, one from 
each of the classes of cubes invariant under integer translations and cube symmetries. For each one 
of these cubes we can construct as before a regular octahedron, obtaining this way a sequence of 
irreducible regular octahedrons. 

2. Some new ingredients and other theoretical facts 

In |12j . we improved and adapted the earlier code for counting all cubes with vertices in 
{0, 1, ..., n} 3 and extended the sequence A098928. In this paper, the usual techniques and ideas are 
going to be the same except some counting procedure that is very efficient in comparison to what 
we had before. We are going to treat this in the general case so, let us suppose that these objects 
can be either equilateral triangles, regular tetrahedrons, cubes or regular octahedrons with vertices 
in Z 3 . For such an object, say O, we can translate it, within Z 3 , to O' that is in the positive octant 
and in such way each plane of coordinates contains at least one vertex of O' . Let us denote by ao 
the number of objects in C m obtained by applying to O' all 48 possible symmetries of the cube 
C m . These symmetries are generated in the following way: first we have symmetries with respect 
to the middle planes and compositions, for example 



(x, y, z) — >• (m — x, y, z), (x, y , z) — ^ (m — x, m — y, z), (x, y,z) — )■ (m — x,m — y,m — z), 

in a total of eight including the identity, then each one of these is coupled with one of the six 
permutations of the variables (Sq). These transformations form a group isomorphic with the group 
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of all 3 x 3 orthogonal matrices having entries ±1 and it is also known as the group of symmetries 
of a cube or of a regular octahedron. It is isomorphic to S4 x Z2. We are going to denote this group 
by S cu be although it is usually known under the name of (extended) octahedral group and denoted 
simply by O h . 

If we think of «o as the cardinality of 

OrbitiO 1 ) := {s(0')\s G S cube }, 

which is, by the first theorem of isomorphism of groups, the same as the cardinality of the group 
factor S cu be/Q, where Q is the subgroup of S cu b e of those symmetries that leave O' invariant. The 
structure of subgroups of S cu b e is known and for each divisor of 48 there is a subgroup of that order. 
Hence, we expect «o to be in the set {1,2,3,4,6,8,12,16,24,48} and most of the time to be 48 
since an arbitrary object O' in C m is unlikely to be invariant under any of the symmetries of S cu b e - 

Then, we denote by a, the cardinality of the set of all the objects counted in «o an d their (all 
possible) integer translations that leave the resulting objects in C m . Also, we denote by f3 the 
objects counted in a which are in {0, 1, ...,m} 2 x {0, 1, ...,m — 1}. Finally, let us denote by 7 the 
number of objects counted in (3 which are in {0, 1, ...,m} x {0, 1, ...,m — l} 2 . Then, we found a 
formula that gives the number of objects obtained from O, under all symmetries and translation 
that leaves the resulting object in {0,1,..., A;} 3 , k > m. 

This fact has been essentially proved in Theorem 2.2 in |10j . The formula that gives this number 
is 

(1) N(0, k) = (k-m + lfa - 3(k - m)(k - m + 1) 2 (3 + 3(fc - m + l)(fc - m) 2 j. 

Let us suppose that the object O can be squeezed within a box of dimensions m x n x p (m > 
n > p), i.e. up to symmetries and translations, O can be transformed to O' fitting snugly into 

B m ,n, P ■= {0, l,...,m} x {0,1,..., n} x {0,1,..., p}. 

We can similarly consider all eight reflections compatible with the box B m)Tl ^ of the form 

(x, y, z) — >■ {m — x, y, z), (x, y,z) — )• (m — x,n — y, z), (x, y,z) — > (m — x,n — y,p — z),etc. 

Let us denote the group of these transformations by Sb- We notice that each one of these 
transformation leaves the object O' inside the box B m , n)V . From case to case, depending of what 
the values m, n and p are, we may have the result of some or all of the permutation transformations 
applied to O' still in B mt n tP . Hence, we will denote by lo(0) the cardinality of the set 
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BoxOrbit(O') := {[ Sl o s 2 ]{0') G £ m , n , P |si G 5 6 ,s 2 G 5 6 }. 

Let us look at an example. Suppose O (equal with O') is the equilateral triangle given by its 
vertices: 

{[0,2, 2], [5, 7,0], [7, 0,1]}. 

We observe that O G -67,7,2- Then one can check that BoxOrbit{0) is the collection of eight 
triangles 

O, {[0,0,1], [2, 7,0], [7, 2, 2]}, {[0,0,1], [2, 7, 2], [7, 2,0]}, {[0,2,0], [5, 7, 2], [7, 0,1]}, 

{[0,5,0], [5, 0,2], [7, 7,1]}, {[0,5, 2], [5, 0,0], [7, 7,1]}, {[0,7,1], [2, 0,0], [7, 5, 2]}, {[0,7,1], [2, 0,2], [7, 5,0]}, 
so uj{0) = 8. It turns out that a {O) = 48, a(O) = 144, p(0) = 40 and 7(C) = 0. Formula © 
becomes 

N(0,k) = 2A(k-l)(k-6) 2 , k>7. 
It turns out the this factorization is not accidental and the following alternative to ([TJ is true. 

Theorem 2.1. Given O, one of the objects mentioned before, and B m ^ n)V the smallest box 
containing a translation of O (m > n > pj, we let u = m — n, v = n — p, and 

A = uj(0){k - m + l)(fc - n + l)(k - p + 1). 

T/ien i/ie number of distinct objects in the cube B^f.^ (k > m), obtained from O by all possible 
integer translations and symmetries is equal to 

A if u = v = 0, 



(2) N(0,k) = { 



3A if u or v is 0, 



6A if u and v > 0. 

Proof. The case u = v = implies w(O) = a (O) = a(O) and /3(C) = 7(C) = because 
there is no room to shift the orbit Orbit(0') inside of -B m , m , m . The formula follows from ([1]). 

Let us look into the case u > and v > 0. We begin by observing that each integer translation 
of the box B mn>p in all possible ways inside B/~ i~k wm gi ve W (C) more copies of O. There is no 
overlap between these copies because neither one of them can be inside of two distinct translations 
of B m7li p. This is due to the minimality of m, n and p. We get A such copies by counting all 
possible translations. Since m, n and p are all distinct, the box -B m ,„, p can be positioned first with 
the biggest of its dimensions along one of the directions given by the axis of coordinates, that is 
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three different ways, and for each such position the next largest dimension can be positioned along 
the two remaining directions. The minimality of m, n and p makes the six different situations 
generate distinct objects. This explains the factor of six that appears in (J2J for this situation. 

In the last case, the box B myTl)P has two of its dimensions the same, so there are only three 
possibilities to arrange the box before one translates it. To see that we get all possible translates 
and symmetries of O by this counting, we can start with one copy O' . Construct the minimum box 
around it. In terms of its position and dimensions, we know in what of the six or three cases we are. 
We transform it into the standard standard position, B m ^ n ^ v , and look at the corresponding object, 
O" . The transformations involved form a group of transformations generated by the permutations 
of the coordinates, the reflections into the axes and integer translations. Every transformation in 
this group, say g = r o a o n with n a permutation, a a reflection or a composition of reflections 
and t a translation, which satisfies g(0) = O" determines a representation (si o S2j{0) = O" with 
s\ G Sb, S2 £ 56 as i n the definition of u>(0). This can be done by taking s 2 = vr and si = tog. This 
is true again because of the minimality of the box B m ^ n ^, i.e. there is only one integer translation 
that takes a reflected box B' mn into -B m , n ,p- B 

This new way of counting is more efficient from a computational point of view because oj is 
simply no bigger than 48, as opposed to the previous situation when a, j3 and 7 could turn out to 
be big numbers and so the number of iterations for computing them would be also large. Roughly 
speaking, this counting factors out fast the problem with the integer translations. 

As an example, let us consider 

OC 2 = {[0,0,1], [0,3, 4], [1,4,0], [3, 0,4], [4, 1,0], [4, 4, 3]}. 

The minimal box here is -84,4,4 and after rotating OC2 in all possible ways (Figure [2] (b)) we get 
uj(OC 2 ) = 4. 






Figure 2 (a): OC2 octahedron 



Figure 2 (b):Four octahedrons in the box 
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The idea of calculations is basically the same as in [12] . in which we have constructed a list of 
irreducible cubes that are used to generate all the other cubes in -Bfc,fe,fc- Here, we are using Theo- 
rem [TTTJ to construct a similar list of irreducible regular hexahedrons. As expected, an irreducible 
regular hexahedron is one whose coordinates cannot be obtained from a strictly smaller one with 
vertices in Z 3 by integer dilations and translations. One simple consequence of Theorem 1 1.1 1 is the 
next corollary. 

Corollary 2.2. The sides of every irreducible regular octahedron in 1? are of the form 

(2k-l)V2 
with k £ N. 

We used the same way of finding the parameterizations of the equilateral triangles as in [12]: 

Q 





Figure 3 (a): Vectors r\ and £ Figure 3 (b): Minimality of parametrization 



(3) 0? = m(-nlf, OQ = n( - (n-m)lf, with (= (d, (1,(2), ~if = (i]i^V2,m), 



(4) 



rac + dbs 



das — bcr 



Xz = r ; 



Vi 



m 



m 



db(s — 3r) + ac(r + s) 
da(s — 3r) — bc{r + s) 



2q 



r + s 



where q = a 2 + b 2 and (r, s) is a suitable solution of 2q = s 2 + 3r 2 that makes all the numbers 



in @ integers. The sides-lengths of AOPQ are equal to d\j2(rrv 2 — ran + n 2 ) and the triangle 
can be completed to a regular tetrahedron (in space) with integer coordinates if and only if k 2 = 
m 2 — mn + n 2 for some k S Z. Related to this fact we have the following proposition. 
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Proposition 2.3. For a prime p > 3, the number of irreducible regular octahedrons in Z 3 (up 
to translations and symmetries) having side lengths equal to py2, is at most ire(p) + 1, where 

,~ ( , A(p) + 24r 2 (3p 2 ) 
(5) 7re(p) = — , 

with A and T defined as in [IT] . 

The exact number of such objects is yet a big mystery to us. 

3. The Maple Code 

We wrote the code using the Maple Software and so we took advantage of the build in commands 
available for a various number of functions. The beginning is pretty standard: 

> restart : with (numtheory) :with(plots) : 

Step 1. The next three procedures calculate all possible values of k, then the parameters m, n, 
and finally the normal vector (a, b, c). The values of k are determined by using a characterization 
theorem for the quadratic form involved here, i.e. all values of k less than n such that k is of the 
form a product of primes of the form 3^ + l,£>2orA; = l. 

> kvalues:=proc(n) 

local i, j ,k,L,a,p,q,r,m,mm; 

L:=0;mm:=floor((n+l)/2) ; 

for i from 2 to mm do 

a:=if actors(2*i-l) ; 

k : =nops (a [2] ) ; r : =0 ; 

for j from 1 to k do 

m:=a[2] [j] [1] ; p:=m mod 3; 
if m=3 then r:=l fi; 
if r=0 and p=2 then r:=l fi; 
od; 
if r=0 then L:=L union {2*i-l};fi; 
od; L:=L union {1}; L:=convert (L,list) ; 
end: 

For example, kvalues(W0) = [1, 7, 13, 19, 31, 37, 43, 7 2 , 61, 67, 73, 79, 7(13), 97]. 

Next, we are interested in the solutions (m,n) of the equation k 2 = m? — mn + n 2 , which are 
primitive in the sense that gcd(m,n) = 1, m > 0, n > 0, and 2m < n. We apply this procedure to 
only those fc's which are the output of kvalues. 
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> listofmn:=proc(k) 
local a,b,i,nx,x,m,n,L; 

x: = [isolve(k"2=m"2-m*n+n"2)] ; 
nx : =nops (x) ; L : ={} ; 
for i from 1 to nx do 

if lhs(x[i] [l])=m then a:=rhs(x[i] [1] ) ; b:=rhs(x[i] [2] ) ; 

else b:=rhs(x[i] [1] ) ; a:=rhs(x[i] [2] ) ; fi; 
if gcd(a,b)=l and a>=0 and b>0 and 2*a<b then L:=L union { [a,b]};f i; 
od;L; 
end: 

A prime of the form 3^ + 1 has a unique primitive decomposition as described, for example, 
listofmn(79) = {[40,91]} since 79 2 = 40 - 40(91) + 91 2 and 79 is a prime. In general, the 
number of solutions is equal to 2 LU ( k > 1 where u>(k) is the number of distinct factors of k which 
are of the form 3£ + 1. This choice of m and n helps identify the irreducible octahedrons: all the 
other solutions (m, n) of the equation k 2 = m 2 — mn + n 2 lead to the same regular octahedron. We 
included the case a = to obtain the output {[0, 1]} for listofmn(l) which is necessary later on. 
In the next procedure, we find only the solutions (a, b, c) of a 2 + b 2 + c 2 = 3d 2 that satisfy 
gcd(a, b,c) = 1 and < a < b < c. The number of such solutions can be calculated precisely 
in terms of the prime factorization of d (see [IT] ) . 

> abcsol:=proc(d) 

local i, j ) k ) m ) u ) x ) y ) sol ) cd;sol:=-[}; 
for i from 1 to d do 

u:=[isolve(3*d~2-i~2=x~2+y~2)] ;k:=nops(u) ; 
for j from 1 to k do 

if rhs(u[j] [l])>=i and rhs(u[j] [2] )>=i then 
cd:=gcd(gcd(i,rhs(u[j] [1] )) ,rhs(u[j] [2])); 
if cd=l then sol:=sol union {sort ( [i,rhs(u[j] [1] ) ,rhs(u[j] [2] )] )};f i;f i; 
od; odjconvert (sol,list) ; end: 

For d = 17, for instance, the procedure abcsol gives the four solutions: [1, 5, 29], [7, 17, 23], [11, 
11, 25] and [13, 13, 23]. We observe that two of them have the property that a = b, in which case, 
we know (see [12]) that the formulae (jl|) simplify quite a bit. 

Step 2. The next seven procedures implement the new way of finding r and s which appear 
in the parametrization (J3J. It is followed by the calculation of the of the parametrization of 
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equilateral triangles ^ and only one regular octahedron is constructed based on Theorem ll.il This 
octahedron is then translated minimally into the positive octant and the minimal cube containing 
it is computed. 

For a prime p of the form 6£ + 1, there exists an unique decomposition p = x 2 + 3y 2 which is 
calculated next. 

> uniquedecomposition:=proc(p) 
local out ,s,outl,out2; 

if p=2 then out := [1,1]; fi; 

if p>2 then s :=isolve(p=x~2+3*y~2) ; 
out 1 : =abs (rhs (s [1] [1] ) ) ; out2 : =abs (rhs (s [1] [2] ) ) ; 
if outl~2+3*out2~2=p then out := [outl,out2] ; else out := [out2,outl] ;f i; 
fi; out; end: 

For p = 2 this procedure has the needed output [1, 1]. As an example, uniquedecomposition(2011) 
[44, 5] since 2011 is a prime and 2011 = 44 2 + 3(5 2 ). The next procedure is calculating the factor- 
ization in Z[\/3i] of a number of the form u + v V3i. 

> f actoroverEisensteinintegers :=proc(u,v) 

local i , N , M , k , a , x , f 1 , f 2 , g , y , y 1 , y2 , L , NN , MM , uu , vv ; 
a : =sqrt (3) *I ; NN : =gcd(u , v) ; uu : =u/NN ; vv : =v/NN ; 
N:=uu~2+3*vv~2;x:=uu+vv*a;M:=if actors (N) ;k:=nops(M[2] ) ; 
for i from 1 to k do 

f 1:=M[2] [i] [1] ; f2:=M[2] [i] [2] ; 
if fl>2 then 

g:=uniquedecomposition(f 1) ; 
else g:=[l,l]; f2:=l; 

fi; 

y :=expand (rationalize (x/(g[l]+a*g [2] ))) ; 
y 1 : =Re (y) ; y2 : =type (yl , integer) ; 

if y2=true then L[i] : = [g[l]+a*g[2] ,f 2] ; else 
L[i] : = [g[l]-a*g[2],f2]; 

fi; 

od; [NN,seq(L[i] ,i=l. .k) , expand (NN*product (L[ii] [l]~L[ii] [2] ,ii=l. .k))] ;end: 
As a simple example here, the following decomposition is obtained for u = 13 and 17: 

(1 + i\/3)(2 + iV3)(5 - 2lV3) = 13 + 171 V3. 
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Perhaps one word of caution is necessary at this point. The decomposition in general is not unique 
in the usual sense since 4 = 2(2) = (1 + v3)(l — \/3). As in the example shown, we go for the 
second representation when something like this happens. We are interested in this decomposition 
because it provides suitable values for r and s that we find next. This turns out to be the greatest 
common divisor between A + Iy/3B and 2q with A = ac, B = bd and q = a 2 + b 2 (see [12] ) . 

> f indgcd:=proc(A,B,q) 

local a,i, j ,f ,f l ) f2 ) common ) m,qq,fac,nf ac,s,rs,P; a:=sqrt (3)*I; 
f :=f actoroverEisensteinintegers(A,B) ;m:=nops(f )-l; 
common : =gcd (f [1] ,2*q) ; qq:=2*q/common~2; P:=common; 
f ac : =if actors (qq) ; nf ac : =nops (f ac [2] ) ; 
for i from 1 to nfac do 

f l:=fac[2] [i] [1] ;f2:=fac[2] [i] [2] ; 
if fl=2 then £2:=l;fi; 

s:=uniquedecomposition(f 1) ; 
for j from 1 to m do 

if s[l]+a*s[2]=f [j] [1] then 

P:=P*(s[l]+a*s[2])~(min(f [j] [2] ,f2)); fi; 
if s[l]-a*s[2]=f [j] [1] then 

P:=P*(s[l]-a*s[2])~(min(f [j] [2] ,f2)); fi; 
od; 
od;P:=expand(P) ;rs:=[Re(P) ,Im(P)/sqrt(3)] ; 
[rs,A*rs[l]+3*B*rs[2] mod 2*q,A*rs [2]-B*rs [1] mod 2*q] ; end: 

In the case we have seen before, where d = 17, if we take a = 1, b = 5 and c = 29, then we have 
A = ac = 29, B = bd = 85 and q = a 2 -\-b 2 = 26. The procedure above gives r = — 1 and s = 7 
and checks that As + 2>Br = (mod 2q) and Ar — Bs = (mod 2q). These two conditions are 
enough to insure that the expressions in @, all give integer values for the coordinates of rj and £. 
In the next procedure we use the r and s determined earlier and construct an irreducible regular 
octahedron in Z 3 , given a vector (a, b, c) and the possible values for (to, n) in the decomposition of 
k 2 = m 2 — mn + n 2 . We are using simple formulae which can be derived easily from Figure 1(b), 
thinking that the point H is the origin, T[l], T[2] and T[3] are the points M, N and L. Then the 
other three vertices are simply given by the addition of every two of these there vectors. Although 
there are six possible equilateral triangles that one may start with in this construction, one can 
see that in the end, essentially the same regular octahedron (up to symmetries and translations) is 
obtained. 
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> f indpar :=proc(a,b,c,mm,nn) 

local q , d , r , s , rs , A , B , k , mx , nx , my , ny , mz , nz , mu , nu , mv , nv , mw , nw , 
u,v,w,x,y,z,T,Rl,R2,DD,E,F; 

q:=a"2+b"2;k:=sqrt(mm"2-mm*nn+nn"2) ;d:=sqrt((a~2+b~2+c~2)/3) ;A:=a*c;B:=b*d; 
rs:=findgcd(A,-B,q) ;r:=rs[l] [2] ; s:=rs[l] [1] ; 
mx : =- (d*b* (3*r+s) +a*c* (r-s) ) / (2*q) ; nx : =- (r*a*c+d*b*s) /q; 

my : = (d*a* (3*r+s) -b*c* (r-s) ) / (2*q) ; ny : =- (r*b*c-d*a*s) /q; mz : = (r-s) /2 ; nz : =r ; 
mu : =nx ; mv : =ny ; mw : =nz ; nu : =nx-mx ; nv : =ny-my ; nw : =nz-mz ; 
u : =mu*mm-nu*nn ; v : =mv*mm-nv*nn ; w : =mw*mm-nw*nn ; 
x : =mx*mm-nx*nn ; y : =my*mm-ny*nn ; z : =mz*mm-nz*nn ; 
Rl : = [ (x+u-2*a*k) /3 , (v+y-2*b*k) /3 , (z+w-2*c*k) /3] ; 
R2 : = [ (x+u+2*a*k) /3 , (v+y+2*b*k) /3 , (z+w+2*c*k) /3] ; 
if Rl[l]=floor(Rl[l]) then T:=[[u,v,w] , [x,y,z] ,R1] ; else 

T:=[[u,v,w] , [x,y,z] ,R2] ; 

fi; 

DD: = [T[1] [1]+T[2] [1] ,T[1] [2] +T [2] [2] ,T[1] [3]+T[2] [3] ] ; 
E: = [T[1] [1]+T[3] [1] ,T[1] [2]+T[3] [2] ,T[1] [3]+T[3] [3]] ; 
F: = [T[2] [1]+T[3] [1] ,T[2] [2]+T[3] [2] ,T[2] [3]+T[3] [3]] ; 
[T[1],T[2],T[3],DD,E,F]; 
end: 

Since for k = 2011 we get listofmn(k) = {[880, 2301]}, we checked to see what regular octahedron 
is obtained for findpar(l, 1, 1,880,2301): [2301, -1421, -880], [880, -2301, 1421], [2401, 100, 1521], 
[3181, -3722, 541], [4702, -1321, 641], and [3281, -2201, 2942]. Since the set a6cso/(2011) has 
336 elements in it, there are essentially at most 337 irreducible regular octahedra in Z 3 with side 
lengths equal to 2011V2 as we pointed out in Proposition 12.31 Next, we have a short function for 
subtracting two vectors U and V. 

subtrv : =proc (U , V) 

> local W; 

W[l] :=U[1]-V[1] ;W[2] :=U[2]-V[2] ;W[3] :=U[3]-V[3] ; [W[l] ,W[2] ,W[3]] ; 
end: 

In order to compare various octahedrons it is easier if they are all translated to the positive octant 
in such a way each plane of coordinates contains at least a vertex of the octahedron. This is 
accomplished with the next routine. 

> tmttopqoctahedron:=proc(T) 
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local i,a,b,c,v,C; 
a:=min(T[l] [1] ,T[2] [1] ,T[3] [1] ,T[4] [1] ,T[5] [1] ,T[6] [1] ) ; 
b:=min(T[l] [2] ,T[2] [2] ,T[3] [2] ,T[4] [2] ,T[5] [2],T[6] [2]); 
c:=min(T[l] [3] ,T[2] [3] ,T[3] [3] ,T[4] [3] ,T[5] [3] ,T[6] [3]); 

v:=[a,b,c] ;C:={subtrv(T[l] ,v) ,subtrv(T[2] ,v) , subtrv(T[3] ,v) , 
subtrv(T[4] ,v) , subtrv(T[5] ,v) , subtrv(T[6] ,v)}; 
end: 

So, for instance, the octahedron mentioned earlier of sides lengths 2011\/2 becomes: [2401, 1521, 
3822], [3822, 2401, 1521], [2301, 0, 1421], [1521, 3822, 2401], [0, 1421, 2301], and [1421, 2301, 0]. 
The smallest cube C m = [0,m] 3 containing an octahedron positioned in the positive octant as 
specified earlier is computed in the following procedure. 

>mscof moctahedron : =proc (Q) 
local a,b,c,T; 
T:=convert(Q,list) ; 
a:=max(T[l] [1] ,T[2] [1] ,T[3] [1] ,T[4] [1] ,T[5] [1] ,T[6] [1]) 
b:=max(T[l] [2] ,T[2] [2] ,T[3] [2] ,T[4] [2] ,T[5] [2] ,T[6] [2]) 
c:=max(T[l] [3] ,T[2] [3] ,T[3] [3] ,T[4] [3] ,T[5] [3] ,T[6] [3]) 
max(a,b,c) ; 
end: 

Step 3. In our construction of these octahedrons we end up with essentially the same octahedron 
if we proceed from a different face of it. To eliminate the possibility of counting an octahedron 
twice or more than one time, we would like to have a way of distinguishing between octahedra 
and so an invariant to translation and symmetries, like the side lengths, will be good. Such an 
invariant is the set of k values which are given by the four pairs of opposite (parallel) faces of a 
regular octahedron. We must have 

I = dik x V2 = d 2 k 2 V2 = d 3 k 3 V2 = d 4 k 4 V2. 

So, knowing £ and the dj's will give us ki's. The set {k 4 ,k2,k 3 ,k 4 } is clearly and invariant to 
translations and the symmetries we have talked about at the beginning of the paper. Hence, 
two octahedra with different sets of k- values will be essentially different. We determine these k- 
values from the first three points given in the procedure findpar. The following calculation finds 
d/gcd(a, b, c), given (a, b, c) such that a 2 + b 2 + c 2 = 3d 2 . 

> unitvector :=proc(U) 
local i,j,k,l,x; 
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i:=U[l];j:=U[2];k:=U[3]; 

l:=gcd(gcd(i, j) ,k) ; x:=(i~2+j~2+k~2)/(3*l~2) ; 
sqrt(x) ;end: 

Then we need a routine to add three vectors. 

> addvec:=proc(U,V,W) 
local X; 

X [1] : =U [1] +V [1] +W [1] ; X [2] : =U [2] +V [2] +W [2] ; X [3] : =U [3] +V [3] +W [3] ; 
[X[l] ) X[2],X[3]];end: 

The next function is multiplying a scalar with a vector. 

> multbyf actorv:=proc(v,k) 

local w; w[l] :=v[l]*k;w[2] :=v[2]*k;w[3] :=v[3]*k; 
[w[l] ,w[2] ,w[3]] ;end: 

The Euclidean distance between to points is needed later. 

> distance:=proc(A,B) 

local C; C:=subtrv(A,B) ; sqrt (C[l] ~2+C[2] "2+C[3] "2) ; end: 

As we said before, in the procedure that follows, T is the list of the first tthree vertices given by 
findpar. 

> f ourkvalues:=proc(T) 
local Nl,N2,N3,N4,x, length; 

length : =distance (T [1] , T [2] ) /sqrt (2) ; 

Nl : =unitvector (addvec (T [1] , T [2] , T [3] ) ) ; 

N2:=unitvector(addvec(T[l] ,T[2] ,multbyfactorv(T[3] ,-3))) 

N3:=unitvector (addvec (T[l] ,T[3] ,multbyfactorv(T[2] ,-3))) 

N4:=unitvector (addvec (T [3] ,T[2] , multbyf actorv(T[l] ,-3))) 
{length/Nl , length/N2 , length/N3 , length/N4} ; 
end: 

For the octahedron in Step 2, we have the set of fc-values, {1, 2011}, as expected since 2011 is a 
prime. 

Step 4. In this step we calculate the orbit of an octahedron T within the reduced cube. The 
procedure orbitboxl is taking care only of the eight possible symmetries. 

> orbitboxl :=proc(T) 

local i,k,Tl ) a ) b ) c ) x ) T2,T3,T4,T5,T6,T7,T8,T9,T10,Tll,T12,T13,T14,T15,T16,T17,T18, 
T19 , T20 , T21 , T22 , T23 , T24 , S , Q ; 
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Q : =convert (T , list) ; 

a:=max(Q[l] [i],Q[2] [1],Q[3] [1],Q[4] [1],Q[5] [1],Q[6] [1] ) 
b : =max (Q [1] [2] , Q [2] [2] , Q [3] [2] , Q [4] [2] , Q [5] [2] , Q [6] [2] ) 
c : =max (Q [1] [3] , Q [2] [3] , Q [3] [3] , Q [4] [3] , Q [5] [3] , Q [6] [3] ) 



Tl 
T2 
T3 
T4 
T5 
T6 
T7 
T8 



=T; 



:={seq([Q[i] [l],Q[i] [2],c-Q[i] [3]] ,i=l. .6)}; 

; ={seq( [Q [i] [1] ,b-Q [i] [2] , Q [i] [3] ] , i=l . . 6) } ; 

: ={seq( [a-Q [i] [1] , Q [i] [2] , Q [i] [3] ] , i=l . . 6) } ; 

;={seq([a-Q[i] [1] ,b-Q [i] [2] ,Q [i] [3] ] , 1=1 . .6)}; 

: ={seq( [a-Q [i] [1] , Q [i] [2] , c-Q [i] [3] ] , 1=1 . . 6) } ; 

: ={seq( [Q [i] [1] ,b-Q [i] [2] , c-Q [i] [3] ] , 1=1 . . 6) } ; 

;={seq([a-Q[i] [1] ,b-Q[i] [2] ,c-Q[i] [3] ] ,1=1 . .6)}; 
S: ={11,12,13, 14,15,16,17,18}; 
S; 
end: 



The next procedure implements Theorem 12.11 in our situation where the objects of interest are 
octahedrons. 



> orbitbox:=proc(T,p) 
local S,Q,TT,a,b,c,m,SS ; 



i,d,u,v,w,y,mm,nn,x,k; 



V 


w} 


;k:=p- 


-d 


i 


=1. 


.6)}; 




i 


=1. 


.6)}; 




i 


=1. 


.6)}; 




i 


=1. 


.6)}; 




i 


=1. 


.6)}; 





=convert CT,list) ; 

=max(Q[l] [1] ,Q[2] [1] ,Q[3] [1] ,Q[4] [1] ,Q[5] [1] ,Q[6] [1]); 
=max(Q[l] [2] ,Q[2] [2] ,Q[3] [2] ,Q[4] [2] ,Q[5] [2] ,Q[6] [2]); 
=max (Q [1] [3] , Q [2] [3] , Q [3] [3] , Q [4] [3] , Q [5] [3] , Q [6] [3] ) ; 
=max ( a , b , c ) ; u : =d-a ; v : =d-b ; w : =d- c ; y : ={u , 
TT[1] :={seq([Q[i] [3] ,Q[1] [2] ,Q[i] [1]] , 
TT[2] :={seq([Q[i] [2] ,Q[i] [3] ,Q[i] [1]] , 
TT[3] :={seq([Q[i] [1] ,Q[i] [3] ,Q[i] [2]] , 
TT[4] :={seq([Q[i] [2] ,Q[i] [1] ,Q[i] [3]] , 
TT[5] :={seq( [Q[i] [3] ,Q[i] [1] ,Q[i] [2]] , 
S:=orbitboxl(T) ; 
for i from 1 to 5 do 

S:=S union orbitboxl (TT [i] ) ; 
od; 
S:=convert(S,list) ;m:=nops(S) ;SS:={}; 
for i from 1 to m do 

If S[l] [1] [l]<=a and S[i] [2] [l]<=a and S[l] [3] [l]<=a and S[i] [4] [l]<=a and S[i] [5] [l]<=a and S[l] [6] [l]<=a 

and S[i] [1] [2] <=b and S[i] [2] [2] <=b and S[i] [3] [2] <=b and S[i] [4] [2] <=b and S[l] [5] [2] <=b and S[i] [6] [2] <=b 
and S[i] [1] [3] <=c and S[i] [2] [3] <=c and S[i] [3] [3] <=c and S[i] [4] [3] <=c and S[i] [5] [3] <=c and S[i] [6] [3] <=c 
and S[i] [1] [1]>=0 and S[i] [2] [1]>=0 and S[i] [3] [1]>=0 and S[i] [4] [1]>=0 and S[i] [5] [1]>=0 and S [i] [6] [1]>=0 
and S[i] [1] [2] >=0 and S[l] [2] [2] >=0 and S[i] [3] [2] >=0 and S [i] [4] [2] >=0 and S[i] [5] [2] >=0 and S[i] [6] [2] >=0 
and S[i] [1] [3] >=0 and S[i] [2] [3] >=0 and S[i] [3] [3] >=0 and S [i] [4] [3] >=0 and S [1] [5] [3] >=0 and S[i] [6] [3] >=0 
then SS:=SS union {S[i]>; fi; 
od; 
mm:=nops(SS) ; 
nn:=(k+u+l)*(k+v+l)*(k+w+l) ; 

if nops(y)=3 then x:=6*mm*nn;f i; 
if nops(y)=2 then x:=3*mm*nn;f i; 
if nops(y)=l then x:=mm*nn;fi; 
[x,SS,nops(SS) , [u,v]] ; end: 



The next two routines calculate the orbit of an octahedron within its minimal cube C m . This 



orbit has at most 48 elements and it is needed in the process of elimination of octahedrons that 
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have already appeared in the construction. In comparison with the previous orbit, it is bigger and 
invariant to all the symmetries. 



> orbitloctahedron:=proc(T) 

local i,k,Tl,a,b,c,x,T2,T3,T4,T5,T6,T7,T8,T9,T10,Tll,T12,T13,T14, 

T15,T16,T17,T18,T19,T20,T21,T22,T23,T24,S,Q,d,al,bl,cl; 

Q:=convert(T,list) ; 

d:=mscofmoctahedron(T) ; 

[QM [2],Q[k] [3],Q[k] [l]],k=1..6)>; 

[Q[k] [l],Q[k] [3],Q[k] [2]],k=1..6)>; 

[Q[k] [l],Q[k] [2],d-Q[k] [3]],k=1..6)>; 

[Q[k] [2],Q[k] [3],d-Q[k] [l]],k=1..6)>; 

[Q[k] [l],Q[k] [3],d-Q[k] [2]],k=1..6)>; 

[Q[k] [l],d-Q[k] [2],Q[k] [3]],k=1..6)>; 

[Q[k] [2],d-Q[k] [3],Q[k] [l]],k=1..6)>; 

[Q[k] [l],d-Q[k] [3],Q[k] [2]],k=1..6)>; 
l([d-Q[k] [l],Q[k] [2],Q[k] [3]],k=1..6)>; 
|([d-q[k] [2],Q[k] [3],Q[k] [l]],k=1..6)>; 
l([d-Q[k] [l],Q[k] [3],Q[k] [2]],k=1..6)>; 
l([Q[k] [l],d-Q[k] [2],d-Q[k] [3]],k=1..6)>; 
l([Q[k] [2],d-Q[k] [3],d-Q[k] [l]],k=1..6)} ; 
l([Q[k] [l],d-Q[k] [3],d-Q[k] [2]],k=1..6)>; 
[([d-q[k] [l],d-Q[k] [2],Q[k] [3]],k=1..6)>; 
l([d-Q[k] [2],d-Q[k] [3],Q[k] [1]] ,k=l. .6)}: 
l ([d-Q[k] [l],d-Q[k] [3],Q[k] [2]],k=1..6)>; 
!([d-Q[k] [l],Q[k] [2],d-Q[k] [3]],k=1..6)>: 
l([d-Q[k] [2],Q[k] [3],d-Q[k] [l]],k=1..6)}: 
l ([d-Q[k] [l],Q[k] [3],d-Q[k] [2]],k=1..6)>; 
l([d-Q[k] [l],d-Q[k] [2],d-Q[k][3]],k=1..6)}; 
|([d-q[k] [2],d-Q[k] [3],d-Q[k] [1]] ,k=l . .6)}; 
l([d-Q[k] [l],d-Q[k] [3],d-Q[k][2]],k=1..6)}; 
S:={T1,T2,T3,T4,T5,T6,T7,T8,T9,T10,T11,T12,T13,T14,T15,T16, 
T17,T18,T19,T20,T21,T22,T23,T24>; 
S ; end : 



Tl 


=T; 


T2 


=-[seq([ 


T3 


={seq([ 


T4 


=-[seq([ 


T5 


=-[seq([ 


T6 


=-[seq([ 


T7 


={seq([ 


T8 


={seq([ 


T9 


=-[seq([ 


Tl( 


):={seq( 


ti: 


. :={seq( 


ti; 


>:={seq( 


ti: 


3:={seq( 


TU 


: :={seq( 


Tl{ 


5:={seq( 


Tit 


5:={seq( 


Tl" 


? :={seq( 


Tl* 


5:={seq( 


Tl< 


):={seq( 


T2( 


):={seq( 


T2: 


. :={seq( 


T2: 


>:={seq( 


T2C 


3:={seq( 


T2' 


: :={seq( 
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> orbitoctahedron:=proc(T) 
local S,Q,T1; 
Q:=convert(T,list) ; 

Tl : ={seq( [Q [k] [3] , Q [k] [2] , Q [k] [1] ] ,k=l . . 6) } ; 
S:=orbitloctahedron(T) union orbit loctahedron (Tl ) ; 
S; end: 

Step 5. At this point we are ready to build a list of minimal, irreducible octahedrons. With 
each entry we keep, in order, d, the value m which is the size of the minimal cube C m , the six 
vertices of the octahedron, the k- values and the corresponding vector (a, 6, c). In the end, another 
list is created in a new procedure, to get the list of all corresponding reducible cubes needed in 
a calculation for a given dimension n. The construction is a little complicated because there are 
octahedrons for which all fe-values are different of 1. 

> ExtendList :=proc(n,N,L: :list , mm , nn , Orb : :array) 
local d,i,sol,nsol,nel,ttpcube,orb: : array, 

LL: :list ,tnel,NL: :list ,cio,0,pO, k, kvalues, m, exception, Int,NN; 
NN:=floor((N+l)/2) ; 
orb:=array(l. .2*NN+1) ; 
nel:=nops(L) ; 
LL:=L; 

k : =sqrt (mm~2-mm*nn+nn~2) ; 
m:=n*k; 
if m<=N then 
sol:=abcsol(n) ;nsol:=nops(sol) ; 
tnel : =nel ; orb : =0rb ; 

for i from 1 to nsol do 
0:=f indpar(sol [i] [1] ,sol[i] [2] , sol [i] [3] ,mm,nn) ;pO:=tmttopqoctahedron(0) ; 
kvalues : =f ourkvalues ( [0 [1] , [2] , [3] ] ) ; 
exception:=evalb(l in kvalues); 
if k=l or exception=f alse then 

Int :=orbitoctahedron(pO) intersect orb[m]; 
cio:=nops(Int) ; 
if cio=0 then 
orb [m] : =orb [m] union orbitbox(pO,N) [2] ; 
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d:=mscofmoctahedron(pD) ; 
tnel:=tnel+l; 

fi; 
fi; 

od; 
LL:=[seq(L[i] ,i=l. .nel) ,seq(NL[j] , j=nel+l. .tnel)] ; 
f i;LL,orb; 
end: 

To generate all the other octahedrons we have to magnify the irreducible ones. 

> multbyf actoroctahedron : =proc (T , k) 
local i,NT,Q;NT:={}; 

Q:=convert(T,list) ; 

for i from 1 to 6 do 

NT:=NT union {multbyf actorv(Q [i] ,k)>; 

od ; NT ; end : 

The procedure ExtendList is used recursively in the next loop to generate the list needed to 
calculate all octahedrons in Cm for a certain value N . 

> ExtendListuptoN:=proc(N) 

local i, j ,k,l,kv,kvn,NN,L,Orb: : array, mn,nmn,E,n,ii; 
kv : =kvalues (N) ; kvn : =nops (kv) ; 
NN:=floor((N+l)/2); 
L:=[] ;0rb:=array(l. .2*NN+1) ; 
for ii from 1 to 2*NN+1 do 

0rb[ii] :={>; 
od; 
for i from 1 to kvn do 
k:=kv[i] ;mn:=listofmn(k) ; nmn : =nops (mn) ; 
for j from 1 to nmn do 
for 1 from 1 to NN+1 do 
n:=2*l-l; 
E:=ExtendList(n,N,L,mn[j] [1] ,mn[j] [2] ,0rb) ; 
L:=E[1]; 
for ii from 1 to 2*NN+1 do 
0rb[ii] :=E[2] [ii] ; 
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od; 
od; 
od; 
od; 
L;end: 

As an example, for n = 20, we obtain for L := ExtendListuptoN (20); 

[[1, 2, {[2, 1, 1], [1, 2, 1], [1, 1, 0], [0, 1, 1], [1, 1, 2], [1, 0, 1]}, {1}, [1, 1, 1]], 
[3, 4, {[4, 1,0], [0,3, 4], [3, 0,4], [1,4,0], [4, 4, 3], [0,0,1]}, {1,3}, [1,1,5]], 
[5, 10, {[8, 5, 7], [4, 10, 4], [1, 5, 8], [4, 0, 4], [7, 5, 0], [0, 5, 1]}, {1}, [1, 5, 7]], 
[7, 12, {[9, 0,4], [12, 8, 9], [0,4, 3], [3, 12, 8], [4, 3, 12], [8, 9,0]}, {1,7}, [1,5, 11]], 
[9, 16, {[3, 9,0], [11, 7, 16], [14, 4, 4], [3, 0,9], [0,12, 12], [11, 16, 7]}, {1,3}, [1, 11, 11]], 
[11, 18, {[7, 18, 13], [11, 0, 1], [18, 7, 13], [15, 15, 0], [0, 11, 1], [3, 3, 14]}, {1}, [1, 1, 19]], 
[13, 24, {[24, 15, 16], [16, 0, 9], [8, 24, 15], [15, 16, 0], [9, 8, 24], [0, 9, 8]}, {1, 13}, [5, 11, 19]], 
[13, 26, {[17, 13, 0], [24, 13, 17], [12, 26, 12], [7, 13, 24], [12, 0, 12], [0, 13, 7]}, {1}, [7, 13, 17]], 
[15, 28, {[15, 28, 13], [0, 19, 1], [20, 12, 0], [20, 9, 21], [5, 0, 9], [0, 16, 22]}, {1, 3}, [5, 11, 23]], 
[17, 24, {[0,3, 4], [0,20, 21], [13, 0,24], [24, 4, 3], [11, 24,0], [24, 21, 20]}, {1}, [1,5,29]], 
[17, 34,{[23, 17, 0], [15, 34, 15], [15, 0, 15], [0, 17, 7], [30, 17, 23], [7, 17, 30]}, {1}, [7, 17, 23]], 
[19, 36, {[0, 12, 12], [34, 24, 24], [11, 17, 36], [11, 36, 17], [23, 0, 19], [23, 19, 0]}, {1}, [5, 23, 23]], 
[19, 30, {[0, 21, 25], [9, 25, 0], [21, 5, 30], [30, 9, 5], [25, 30, 21], [5, 0, 9]}, {1, 19}, [1, 11, 31]]] 
Let us observe that although there are four primitive solutions for 

abcsol(19) = [[5, 23, 23], [11, 11, 29], [13, 17, 25], [1, 11, 31]] 

we get only two essentially different octahedrons of side lengths 19\/2. We also need to take into 
account the reducible octahedrons. 

> L:=ExtendListuptoN(100) : 
ExtendListuptoNmultiples:=proc(N,L) 
local i, j ,x,lc,m,mm,d,dd,C,CC,LL; 
m : =nops (L) ; i : =1 ; LL : ={} ; 
while i<=m do 
d:=L[i] [2]; 

if d<=N then 

mm:=floor(N/d) ;C:=L[i] [3] ;j:=2; 
while j<=mm do 

CC : =multbyf act or octahedron (C , j ) ; 
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dd : =d* j ; lc : =nops (LL) ; 

LL:=LL union { [L[i] [1] *j ,dd,CC,L[i] [4]] >; 

od; 

fi; 
i:=i+l; 
od; 
convert (LL, list) ; 
end: 

> LL:=ExtendListuptoNmultiples(20,L) : 

The program is taking the previous list and inflates only the octahedrons needed: 

LL := [[3, 6, {[3, 3, 6], [3, 3, 0], [3, 6, 3], [0, 3, 3], [6, 3, 3], [3, 0, 3]}, {1}], 

[2, 4, {[2, 2, 4], [2, 2,0], [2, 4, 2], [0,2, 2], [4, 2, 2], [2,0,2]}, {1}], 

[7, 14, {[7, 7, 14], [7, 7, 0], [7, 14, 7], [0, 7, 7], [14, 7, 7], [7, 0, 7]}, {1}], 

[6, 12, {[6, 6, 12], [6, 6, 0], [6, 12, 6], [0, 6, 6], [12, 6, 6], [6, 0, 6]}, {1}], 

[5, 10, {[5, 10, 5], [0, 5, 5], [5, 5, 10], [5, 5, 0], [10, 5, 5], [5, 0, 5]}, {1}], 

[4, 8, {[4, 0,4], [4, 4,0], [4, 4, 8], [4, 8, 4], [0,4, 4], [8,4, 4]},{1}], 

[10, 20, {[10, 10, 20], [10, 10, 0], [10, 20, 10], [0, 10, 10], [10, 0, 10], [20, 10, 10]}, {1}], 

[9, 18, {[9, 9, 18], [9, 9, 0], [9, 18, 9], [0, 9, 9], [18, 9, 9], [9, 0, 9]}, {1}], 

[8, 16, {[8, 8, 16], [8, 8, 0], [8, 16, 8], [0, 8, 8], [16, 8, 8], [8, 0, 8]}, {1}], 

[15, 20, {[0, 15, 20], [0, 0, 5], [15, 0, 20], [5, 20, 0], [20, 20, 15], [20, 5, 0]}, {1, 3}], 

[12, 16, {[0, 12, 16], [4, 16, 0], [16, 4, 0], [0, 0, 4], [12, 0, 16], [16, 16, 12]}, {1, 3}], 

[9, 12, {[0, 9, 12], [0, 0, 3], [9, 0, 12], [12, 12, 9], [3, 12, 0], [12, 3, 0]}, {1, 3}], 

[6, 8, {[2, 8,0], [8, 8, 6], [8, 2,0], [0,6, 8], [6, 0,8], [0,0, 2]}, {1,3}], 

[10, 20, {[8, 0, 8], [14, 10, 0], [8, 20, 8], [0, 10, 2], [16, 10, 14], [2, 10, 16]}, {1}]] 

Step 6. Finally we are adding up the contribution of each octahedron located in the union of 
the two lists created earlier. 

> addupnew:=proc(N,L,LL) 

local i, j ,k,nc,mm,m,d,C,CC,x,dd,nt; nc:=0; m:=nops(L); i:=l; 
while i<=m do 
d:=L[i] [2]; 

if d<=N then 
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x:=orbitbox(L[i] [3] ,N) [1] ; 
nc:=nc+x; 
fi; i:=i+l; 
od; 
m:=nops(LL) ; 
i:=l; 

while i<=m do 
d:=LL[i] [2]; 
if d<=N then 

x:=orbitbox(LL[i] [3] ,N) [1] ; 
nc:=nc+x; 
fi; i:=i+l; 
od; 
nc;end: 

The next command produces the sequence we are looking for. 
>NO:=[seq([k,addupnew(k,L,LL)] ,k=l. .100)] ; 

So, the first one hundred terms of A178797 are: 

NO := [[1, 0], [2, 1], [3, 8], [4, 32], [5, 104], [6, 261], [7, 544], [8, 1000], [9, 1696], [10, 2759], 
[11, 4296], [12, 6434], [13, 9352], [14, 13243], [15, 18304], [16, 24774], [17, 32960], [18, 43223], [19, 
55976], [20, 71752], [21, 90936], [22, 113973], [23, 141312], [24, 173436], [25, 210960], [26, 254587], 
[27, 305000], [28, 364406], [29, 432824], [30, 511421], [31, 600992], [32, 702556], [33, 817200], [34, 
946131], [35, 1090392], [36, 1251238], [37, 1430072], [38, 1629391], [39, 1850064], [40, 2094276], [41, 
2363616], [42, 2659813], [43, 2984600], [44, 3341660], [45, 3731720], [46, 4156689], [47, 4618480], 
[48, 5119292], [49, 5661600], [50, 6248705], [51, 6882808], [52, 7568126], [53, 8306520], [54, 9104339], 
[55, 9962320], [56, 10888762], [57, 11882896], [58, 12949661], [59, 14090952], [60, 15311286], [61, 
16613736], [62, 18001975], [63, 19479680], [64, 21052826], [65, 22724576], [66, 24500175], [67, 
26383240], [68, 28387456], [69, 30510616], [70, 32758963], [71, 35136544], [72, 37656214], [73, 
40317328], [74, 43125329], [75, 46085496], [76, 49207224], [77, 52493112], [78, 55954267], [79, 
59592272], [80, 63415296], [81, 67428832], [82, 71642127], [83, 76059704], [84, 80701546], [85, 
85565064], [86, 90662451], [87, 95997360], [88, 101592122], [89, 107443264], [90, 113561009], [91, 
119951832], [92, 126644136], [93, 133629672], [94, 140916757], [95, 148513712], [96, 156444624], [97, 
164706400], [98, 173308509], [99, 182260568], [100, 191575248]] 
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